Superheating field of superconductors within Ginzburg-Landau theory 
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We study the superheating field of a bulk superconductor within Ginzburg-Landau theory, which 
is valid near the critical temperature. We calculate, as functions of the Ginzburg-Landau parameter 
k, the superheating field H s ^ and the critical momentum k c characterizing the wavelength of the 
instability of the Meissner state to flux penetration. By mapping the two-dimensional linear stability 
theory into a one-dimensional eigenfunction problem for an ordinary differential equation, we solve 
the problem numerically. We demonstrate agreement between the numerics and analytics, and 
show convergence to the known results at both small and large k. We discuss the implications of 
the results for superconducting RF cavities used in particle accelerators. 

PACS numbers: 74.25.Op 



I. INTRODUCTION 



One of the primary features of superconductivity is 
the Meissner effect — the expulsion of a weak magnetic 
field from a bulk superconducting material [l| . For suffi- 
ciently large magnetic fields, the Meissner state becomes 
unstable, and the system undergoes a phase transition. 
The exact nature of the transition depends on the so- 
called Ginzburg-Landau parameter, k = A/£, where A is 
the London penetration depth and £ the superconducting 
coherence length. Type I superconductors, characterized 
by small k, transition from the Meissner state into a nor- 
mal metal state for magnetic fields above the thermody- 
namic critical field, H c . Type II superconductors, with 
larger k, instead transition into a superconducting state 
with vortices above the first critical field H c \ . This state 
is stable up to a second critical field H C 2 , above which the 
metal becomes normal. For any superconductor, how- 
ever, the Meissner superconducting state is metastable, 
persisting up to the superheating field , well above H c 
or H c i (for type I and II superconductors, respectively). 
The main goal of the present work is the calculation of 
-ff s h as function of k for superconductors near the critical 
temperature T c where Ginzburg-Landau theory is appli- 
cable (we remind that within Ginzburg-Landau theory, 
the transition from type I to type II superconductors is 
at k = l/y/2). 

The metastability of the Meissner state is of interest 
in the design of resonance RF cavities in particle accel- 
erators, where H s h places a fundamental limit on the 
maximum accelerating field Q. As type II supercon- 
ducting materials are being considered in cavity designs, 
a precise calculation of in this regime is of value. 
One must note, however, that operating temperatures 
of superconducting RF cavities are well below the criti- 
cal temperature T c and that at these low temperatures 
Ginzburg-Landau theory is not quantitatively valid. The 
numerical techniques developed here are also being used 
within the Eilcnbcrgcr formalism to address these lower 
temperatures [3( . Using this formalism, the limit k — > 00 
was studied in Ref. [4] for arbitrary temperature. 



Much work has already been done in calculating the su- 
perheating field within Ginzburg-Landau theory [514121]. 
The problem is formulated as follows: the superconduc- 
tor occupies a half space with a magnetic field applied 
parallel to the surface. The order parameter and vector 
potential are functions of the distance from the surface 
and can be found by solving a boundary value problem of 
ordinary differential equations. The superheating field is 
then the largest magnetic field for which the correspond- 
ing solution is a local minimum of the free energy. For 
small values of ft, the superheating field corresponds to 
the largest magnetic field for which a nontrivial solution 
to the Ginzburg-Landau equations exist [l2[ , as the insta- 
bility does not break translational invariance. However, 
as k increases the one-dimensional solution is unstable to 
two-dimensional perturbations, resulting in a lower esti- 
mate of H s h as first shown in Ref. Q ■ The task at hand 
is to find which perturbations destroy the Meissner state 
and at which value of the applied magnetic field they first 
become unstable. 

The calculation of iJ s h is therefore a linear stability 
analysis of the coupled system of superconducting order 
parameter and vector potential. For a given configura- 
tion, we study its stability to arbitrary two-dimensional 
perturbations by considering the second variation of the 
free energy: if the second variation is positive defi- 
nite for all possible perturbations then the solution is 
(meta)stable. The second variation can be expressed as 
a Hermitian operator acting on the perturbations, so it 
is sufficient to show that the eigenvalues of this oper- 
ator are all positive. By expanding the perturbations 
in Fourier modes parallel to the surface, the eigenvalue 
problem can once again be translated into a boundary 
value problem of an ordinary differential equation. The 
eigenvalues now depend upon the wave-number of the 
Fourier mode, but can otherwise be solved in the same 
way as the Ginzburg-Landau equations. The superheat- 
ing field is then the largest applied magnetic field for 
which the smallest eigenvalue is positive for all Fourier 
modes. 

The present stability analysis is more challenging than 
many such calculations, as the instability destabilizes an 
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interface with a pre-existing depth-dependence of field 
and superconducting order parameter. As described 
above, we map the partial differential equation for the 
unstable mode into an eigenvalue analysis for a family of 
one-dimensional ordinary differential equations (as orig- 
inally suggested, but not implemented, in Ref. [5]). This 
technique could be useful in a variety of other linear sta- 
bility calculations (T3l - fl5l |. replacing thin interface ap- 
proximations with a microscopic depth-dependent treat- 
ment of the destabilizing interface. 

The paper is organized as follows: in the next section 
we present the Ginzburg-Landau free energy and the dif- 
ferential equations to be studied for the stability analysis. 
In Sec. IIIII we give some details about the numerical cal- 
culations and our main results. In Sec. lIVI we discuss the 
implications of the results for accelerator cavity design 
and outline future research directions. In Appendices we 
derive analytic formulas, valid at large k, which we com- 
pare against the numerics. 

II. GINZBURG-LANDAU THEORY AND 
STABILITY ANALYSIS 

The Ginzburg-Landau free energy for a superconductor 
occupying the half space x > in terms of the magnitude 
of the superconducting order parameter / and the gauge- 
invariant vector potential q is given by 

Jx>a 1 1 

+/V + (H a -AVxq) 2 }, (1) 

where H Q is the applied magnetic field (in units of 
y/2H c ), £ is the Ginzburg-Landau coherence length, and 
A is the penetration depth. Note that after choosing 
the unit of length, the only remaining free parameter in 
the theory is the ratio of these two characteristic length 
scales, the Ginzburg-Landau parameter n — A/£. The 
magnetic field inside of the superconductor is given by 
H = AV x q. 

We take the applied field to be oriented along the z- 
axis H a = (0,0, H a ), and the order parameter / = f(x) 
to depend only on the distance from the superconductor's 
surface. We have assumed that the order parameter is 
real and further parametrize the vector potential as q = 



(0, q(x), 0), which fixes the gauge. The Ginzburg-Landau 
equations that extremize T with respect to / and q are 

er q 2 f +.f- f = o, 

AY - fq = 0, U 

and with our choices H = Xq' '. Hereafter we use primes 
to denote derivatives with respect to x. 

The boundary conditions at the surface derive from 
the requirement that the magnetic field be continuous, 
q'(0) = H a /X, and that no current passes through the 
boundary, /'(0) = 0. We also require that infinitely 
far from the surface the sample is completely super- 
conducting with no magnetic field, giving us f(x) — > 1 
and q(x) — > as x — > oo. In the limits k — ¥ and 
k — »• oo, Eqs. (HJ can be explicitly solved perturbatively, 
see Ref. [12j ] and Appendix [3 respectively. For arbitrary 
k they can be solved numerically via a relaxation method, 
as we discuss in Sec. IIIII 

For a given solution (/, q) we consider the second varia- 
tion of T associated with small perturbations / — > / + Sf 
and q — » q + Sq given by 

S 2 T = f ^£ 2 (VS/) 2 + 4/5/q-5q + / 2 <5q 2 

Jx>0 L 

(3/ 2 +q 2 -lM/ 2 + A 2 (Vxiq) 2 }. (3) 

If the expression in Eq. ((3]) is positive for all possible per- 
turbations, then the solution is stable. Since our solution 
(/, <5q) depends only on the distance from the boundary 
(and is therefore translationally invariant along the y and 
z directions), we can expand the perturbation in Fourier 
modes parallel to the surface. As shown in Ref. we 
can restrict our attention to perturbations independent 
of z and write 

Sf(x,y) = S f(x) cos ky, 

Sq(x, y) = (Sq x sin ky, Sq y cos ky, 0), 

where k is the wave-number of the Fourier mode. The 
remaining Fourier components (corresponding to replac- 
ing cos — > sin and vice- versa in Eq. |4]) are redundent as 
they decouple from those given in Eq. 2] and satisfy the 
same differential equations derived below. 

After substituting into the expression ([3]) for the sec- 
ond variation and integrating by parts, we arrive at 



8 2 F- 



dx ( Sf Sq y Sq x ) 



2 d z 
tlx 2 



+ 3/ 2 

Vq 




i 2 k 2 -l 2fq 



X 2 k 



dx 



f 2 + X 2 k 2 




(5) 



r 



The matrix operator in Eq. (J5j> is self-adjoint, and the second variation will be positive definite if its eigenvalues 
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are all positive. In the eigenvalue equations for this op- 
erator, the function 5q x can be solved for algebraically. 
The resulting differential equations for Sf and 8q y are 



fSf" + (3/ 2 + q 2 - 1 + ek 2 )Sf + 2fqSq y 



ESf, 
(6) 



and 



dx 



P-E 



P + X 2 k 2 -E 



f 2 S% 



2fqSf = ESq 



(7) 

where E is the stability eigenvalue. Note that by decom- 
posing in Fourier modes, we have transformed the two- 
dimensional problem into a one-dimensional eigenvalue 
problem. Numerically, it can be solved by the same re- 
laxation method as the Ginzburg-Landau equations - see 
Sec. Mil The boundary conditions associated with the 
eigenvalue equations derive from the same physical re- 
quirements previously discussed: we require 8f'(0) = 0, 
since no current may flow through the boundary, and 
Sq' y (0) — 0, since the magnetic field must remain contin- 
uous. Additionally, we require 5f(x) and Sq y (x) 
as x — ¥ oo. There is also an arbitrary overall normaliza- 
tion, which we fix by requiring 5f(0) = 1. 

The stability eigenvalue will depend on the solution of 
the Ginzburg-Landau equations, i.e., the applied mag- 
netic field H a , and the Fourier mode k under considera- 
tion. The problem at hand is to find the applied magnetic 
field and Fourier mode for which the smallest eigenvalue 
first becomes negative, which is the case if the following 
two conditions hold: 



E = 0, 



dE 
~dk 



0. 



(8) 



The value of the magnetic field at which these conditions 
are met is the superheating field H s h , and the correspond- 
ing wave-number is known as the critical momentum k c . 
In the next section we discuss in more detail the numer- 
ical approach used to calculate these two quantities. 



III. NUMERICAL RESULTS 

As explained in the previous section, the calculation 
of the superheating field comprises two main steps: (1) 
solving the Ginzburg-Landau equations and (2) solv- 
ing the eigenvalue problem ©-(0) with conditions (JSJ). 
To solve these equations we employ a relaxation method. 
The basic scheme is to replace the ordinary differential 
equations with a set of finite difference equations on a 
grid. From an initial guess to the solution, the method 
iterates using Newton's method to relax to the true so- 
lution [lj|. The grid is chosen with a high density of 
points near the boundary, with the density diminishing 
approximately as the inverse distance from the bound- 
ary. This is similar to the scheme used by Dolgert et al. 
[12j ] to solve the Ginzburg-Landau equations for type I 
superconductors. 



For k near the type I/II transition, the relaxation 
method typically converges without much difficulty. In 
the limiting cases that k becomes either very large or 
small, however, the grid spacing must be chosen with 
care to achieve convergence. The eigenfunction equa- 
tions are particularly sensitive to the grid choice. This 
is not surprising, since in either limit there are two well- 
separated length scales. For example, using units A = 1 
and £ = 1/k, we find that a grid with density 



p{x) = 



150k 
1 + 25kx 



(9) 



leads to convergence for k as high as 250. The grid points 
are then be generated recursively x*+i = Xf+1/ p(xi) with 
xq = 0. We find that if the grid is not sufficiently sparse 
at large x, the relaxation method fails, presumably due to 
rouding errors. On the other hand if it is too sparse, the 
finite difference equations poorly approximate the true 
differential equation. Fortunately, the method converges 
quickly, allowing us to explore the density by trial and 
error, as we have done to get Eq. [§] 

In solving Eq. [3J if a sufficiently large value for the ap- 
plied magnetic field is used, there may not be a nonzero 
solution to the Ginzburg-Landau equations and the re- 
laxation method will often fail to converge, indicating 
that the proposed H a is above the actual superheating 
field. In practice, therefore, it is more convenient to re- 
place the boundary condition q'(0) = XH a with a con- 
dition on the value of the order parameter /(0) = A, 
which then implicitly defines the applied magnetic field 
as a function of A, H a (A). This has the advantage that 
H a (A) is a differentiable function of A, as is the stabil- 
ity eigenvalue, improving the speed and accuracy of the 
search for the superheating field. The drawback to this 
approach is that H a (A) is not single- valued, with an un- 
stable branch of solutions as illustrated in Fig.[TJ For the 
problem at hand, this turns out to be straightforward to 
address since we determine the stability of each solution 
in the second step. To achieve the conditions in Eq. (JSJ), 
we vary both the Fourier mode, k, and the value of the 
order parameter at the surface A. 

The results of the procedure described above are sum- 
marized in Figs. [2][5l where we also compare them with 
analytical estimates which, for large values of k, are 
derived in Appendices. In Fig. [5] we plot the numer- 
ically calculated superheating field as a function of k 
(solid line). The vertical line at k c ~ 1.1495 separates 
the regimes of one-dimensional (ID, k = 0) and two- 
dimensional (2D, k ^ 0) critical perturbations. We have 
checked the value of k c both by assuming 2D pertur- 
bations and finding when their critical momentum goes 
to zero and by assuming ID perturbation and finding 
when the coefficient of the term quadratic in momen- 
tum in the second variation of the free energy vanishes; 
these methods lead to the same value within our numer- 
ical accuracy of 10~ 4 . Our value of k c is higher than 
previous estimates, which ranged from 0.5 @ to 1.10 @ 
and 1.13(±0.05) [10(. k c is larger than the boundary 
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H a /\/2H c 

FIG. 1: Solving for the 1-D Ginzburg-Landau super- 
heating field. By fixing value of the order parameter at the 
surface, we implicitly define the applied magnetic field H a {A) 
at the surface. This definition produces a branch of unsta- 
ble solutions, but guarantees that our equations will have a 
solution for all guesses of A. The "nose" of the curve occurs 
at the , the superheating field ignoring two-dimensional 
fluctuations, and is the largest H a for which a nontrivial so- 
lution to Eq. ([2]) can be found. This example was calculated 
for k = 1, and « 0.9. 
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FIG. 3: Comparison between the numerical and asymptotic 
critical momentum k c , Eq. (|13p . The approximate behavior 
of fc c near re c (dot-dashed) is given in Eq. (|12p . 

proximation to the actual H s h, with deviation of less than 
about 1.5 %. In the opposite case k > k c , 2D perturba- 
tions are the cause of instability and the superheating 
field is approximately given by [lOJ (dashed line) 

H sh (n) VTO 0.3852 

V2H C ~ 6 + ■ [ ' 

Equation (fTTj) . derived in Appendix [B] is also a good 
approximation, deviating at most about 1 % from the 
numerics. Therefore, our numerics show that the simple 
analytical formulas for in Eqs. (fT0|) and (jlll) can 
be used to accurately estimate the superheating field for 
arbitrary value of the Ginzburg-Landau parameter, when 
used in their respective validity regions. 

In Fig. [3] we show the numerical result for the critical 
momentum k c versus k. We see that k c — ¥ as k — > k c 
from above. Near k c , the behavior of k c is reminiscent 
of that of an order parameter near a second-order phase 
transition: 



FIG. 2: Numerically calculated H s h and corresponding 
analytical approximations [Eqs. (|10[) and (| 11 f)] versus the 
Ginzburg-Landau parameter re. 

k = 1/ v2 separating type I from type II superconduc- 
tors. Type II superconductors for which k < n c become 
unstable via a spatially uniform invasion of magnetic flux. 
Additionally, we find that superheated type II supercon- 
ductors with k < 0.9192 can transition directly into the 
normal state since the corresponding i? s h is larger than 
the second critical field H C 2. 

For k < k c the instability is due to ID perturbations. 
In this regime, the Pade approximant 

H sh (n) _ 2 _ 3/4 _ 1/2 1 + 4.6825120k + 3. 3478315k 2 
y/2H c ~ K 1 + 4.0195994k + 1.0005712k 2 

_ (10) 
derived in Ref. [12[ (dot-dashed line) gives a good ap- 



k c ~ 1.2-y/K- k c , (12) 

where the prefactor has been estimated by fitting the 
numerics. The dashed line is the asymptotic formula [Toj 
(see also Appendix |B|) 

Xk c w 0.9558k 3/4 (13) 

which captures correctly the large-K behavior. 

In Fig. |4] we present a typical solution for the order 
parameter near the surface at H = H s h for a large value 
of k, along with the analytic approximations presented 
in Appendix [A} The zeroth order approximation /o, 
Eq. (|A2[) . fails near the surface, as it does not satisfy 
the boundary condition /'(0) = 0. On the other hand, 
including the first order correction in 1/k, Eq. (|A15[) . 
leads to excellent agreement with the numerics. Finally, 
in Fig. [S]we show a typical example of the depth depen- 
dence of the perturbation Sq y at the critical point where 



5 



0.80r 

0.79 
0.78- 
0.77- 
0.76- 

°-4- 



— numerical 

fo 

— - /o + /i/k 



2 3 



10' 



a/A 
2tt/AAv. 



K 



W 2 



FIG. 4: Profile of the order parameter at for k = 50 
together with analytic approximations given in Eqs. (|A2jl and 



FIG. 6: The wavelength of the critical perturbation (2n/k c ) 
and the Abrikosov vortex spacing (a) calculated at the su- 
perheating field both vanish at large n, although the former 
diminishes much more quickly. 
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x/X 

FIG. 5: Numerical profile (solid line) of the critical perturba- 
tion Sq y determining H s h for ft — 50 compared to the large-K 
perturbative result (dashed), Eq. (|B22[1 . 



the solution first becomes unstable. We find again good 
agreement between numerical and perturbative calcula- 
tions. 

It is interesting to compare the wavelength of the crit- 
ical perturbation, 27r/fc c , with the Abrikosov spacing a 
for the arrangement of vortices at the superheating field: 
the naive expectation is that the initial flux penetrations 
represent nuclei for the final vortices. Kramer argues that 
this picture is incorrect since the initial flux penetrations 
do not have supercurrent singularities and do not carry 
a fluxoid quantum [j| . 

We find the numerical discrepancies between the two 
lengths further support Kramer's argument. In the 
weakly type-II regime (k ~ 1), the initial flux penetra- 
tion is from infinitely long wavelengths (k c = 0); in con- 
trast, the final vortex state has a very high density, since 
H s h ~ H c2 (l7j |. In the strongly type-II limit (k — > oo) 
both the inverse momentum and the Abrikosov spacing 
(evaluated at the superheating field) vanish, but at differ- 
ent rates, with l/k c ~ k -3 / 4 while the Abrikosov spacing 
see Fig. [5] These results suggest 



that there is no immediate connection between the ini- 
tial penetration and the final vortex array. A dynamical 
simulation could explore the transition between the ini- 
tial penetration and the final vortex state, similar to that 
done by Frahm et. al for the transition from the normal 
state to the vortex state [l9| . 



IV. SUMMARY AND OUTLOOK 

In this paper, we have numerically calculated within 
Ginzburg-Landau theory the superheating field H s h of su- 
perconductors by mapping the linear stability threshold 
onto an eigenfunction problem, and we have shown that 
analytic approximations are in good agreement with the 
numerical results. The technique of mapping the linear 
stability problem onto a one-dimensional eigenfunction 
problem is potentially a useful technique, and we hope 
others find useful applications of the methods described 
here. 

One of the primary motivations for this work is the ap- 
plication to RF cavities in particle accelerators, where the 
maximum accelerating field is limited by H s h . While the 
results presented here provide good estimates of H s h for 
many materials of interest near the critical temperature 
T c , we emphasize that the operating temperature of these 
cavities are typically well below T c , where Ginzburg- 
Landau theory is not quantitatively accurate. The tech- 
niques presented here can be applied to Eilenberger the- 
ory to more accurately determine i? s h at low tempera- 
tures. The Eilenberger approach has already been used 
[3| to evaluate H s h (T) at any temperature in the infinite 
k limit for clean superconductors, and work is in progress 
to address low temperatures for finite k Q. 
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Appendix A: Order parameter and vector potential 
in the large k limit 

In this Appendix, we derive solutions to the Ginzburg- 
Landau equations Eq. @ valid in the large-K limit. For 
convenience, we work in units A = 1, £ = l/n. As a first 
step, we consider the limit k — > oo. Then Eqs. [2] reduce 
to 



3o = foQo, 
= /o (/o 2 - 1 



(Al) 



with solution 



?o0) 



V2 



cosh(a; + £) ' 
fo(x) = yjl-q*(x), 



(A2) 



where the parameter t is determined by the field at the 
surface via 



9o(0) 



V2sinh£ 
cosh 2 1 



(A3) 



The above solution satisfies the boundary conditions 
at infinity, but it cannot satisfy the boundary condition 
for / at the surface. An approximate solution, valid at 
finite but large k, which satisfies all boundary conditions 
can be obtained by boundary layer theory. We follow the 
approach of Ref. [ljj , so we only sketch the steps of the 
calculation. Note that away from the thermodynamic 
critical field, the scaling is different than that used in 
Ref. [ll|: there the expansion is in powers of n~ a and 
the inner variable is X — n a x with a = 2/3, here we use 
a = 1: 



<7o 



1 

— < 

K 
1 



-qi 



(A4) 



/ = So + -h + 



Substituting into Eqs. ©, we find the following "outer 
layer" equations for q± and f\. 



^ = 2/o/igo + ZoV 
= /i(3/ 2 -l + g 2 )+2/o<Zo<Zi 



(A5) 



which have the simple solutions fi = qi =0. For the 
inner layer, we introduce the variable X = kx and find 



and 



fo =/o(/o 2 -l + 9o 2 

q'o = o, 



fi = h (3/o 2 - 1 + qf) + 2/oMi, 
<7i = 0, 



(A6) 



(A7) 



where we use tildes to denote functions of the inner vari- 
able X. Equations IA6I have constant solutions 



'in 



= -b, f = v 7 ! - & 2 , 



(A8) 



while from the second of Eqs. (|A7j) and the boundary 
conditions we get 



& = H a X . 
Then the first of Eqs. (|A7[) becomes 



(A9) 



/{' = 2(1 - b 2 )h - 2bVl - bm a X, (A10) 
with solution 



h = 



_bHa_ 



x + Ae -V2VT^x + Be V2VT^x (An) 



with A, B integration constants. Since / tends to a con- 
stant far from the surface, we set B = 0. Vanishing of 
the derivative at the surface then fixes 



.4 



bH n 



V2(l - & 2 ) 



(A12) 



Next, we match the inner and outer solutions. Com- 
paring Eqs. (|X2|) and (|A"8|) we get 



V2 
cosh< 



(A13) 



We can express b in terms of the applied field using 
Eq. ££21) to find 



l-v/l-2ff2. 



(A14) 



Then, since fi=qi= 0, we need to compare the linear 
order expansion of Eqs. (|A2j) at small x with qi/n and 
/i/k at large X = kx. Using Eqs. (|A3j) and (| A13|) . we 
find that the inner and outer solutions match. Finally, 
the uniform approximate solution is 



q(x) = qo(x), 



1 bH a a _^^ KX (A15) 



V2(l-fo 2 ) 

with corrections of order 1/k 2 . In Fig. 0] we compare the 
second of Eq. (|A15|) to numerics. 
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Appendix B: Superheating field in the large-K limit 

The calculation of the superheating field i? s h as a func- 
tion of k for stability with respect to one-dimensional 
perturbations (i.e., k = 0) can be found in Ref. [Hj for 
k — y and Ref. [llj for « — > oo. The latter calculation, 
however, is of little physical relevance, as the actual in- 
stability at sufficiently large k is due to two-dimensional 
perturbations. Here we present for completeness (albeit 
in a different form) Christiansen's perturbative calcula- 
tion [lfj of the true superheating field iJ s h(«) for k 3> 1. 

Our starting point is the following expression for the 
"critical" second variation of the thermodynamic poten- 
tial as functional of perturbations Sf, 5q y , and momen- 
tum k [see also Eq. (10) in Ref. 



S 2 T- 



{ [3/ 2 + q 2 - 1 + {k/n) 2 } Sf 2 + K - 2 Sf' 2 

ifqSfSq y + f 2 Sql + (f 2 + k^fSqty 

(Bl) 



dx 



It is straightforward to check that variation of this func- 
tional with respect to / and q y leads to Eqs. ©-([7]) with 
E = and rescaled units A = 1. Kramer estimated that 
the critical momentum k oc yOt. While we will show 
that this is not the correct scaling, this form suggests to 
rescale lengths by l/v^ by defining x — w/ ^fh~: 



S 2 T = 



dw 



{[3f 2 + q 2 -l + (k/n) 2 ] Sp + n-'Sf 2 

+ 4fqSfSq v + f 2 5q 2 y + (f + k 2 )' 1 } 2 nSq' y 2 ) 

(B2) 



where now prime is derivative with respect to w. (Note 
that although k has units of inverse length, it is momen- 
tum parallel to the surface, and therefore does not scale 
with x.) 

Minimization with respect to k leads to the equation 



k / dw 



Sf 2 



nf 2 



M 2 



if 2 



(B3) 



Assuming k 3> 1, we can neglect f 2 < 1 in the denomi- 
nator and find 



fc 4 / dw Sf 2 = k 3 dw f 2 5q' y ' 



(B4) 



which shows that (if our length rescaling is correct) the 
proper scaling for the critical momentum is k oc k 3 / 4 . 
If this is true, then {k/n) 2 oc an d n/k 2 oc X/^JH, 

which shows that the next to leading order terms in curly 
brackets in Eq. (|B2[) are proportional to l/^fn. There- 
fore, terms of order 1/k can be neglected and, in par- 
ticular, we can neglect k Sf and use everywhere the 
lowest order solution for / and q, Eq. (|A2I) . Hence the 



approximate functional in the large- n limit is 

r d j [2f Q + ( k / K ) 2 ] sf 2 

VkJo I. ( B5 ) 
+ 4f q Q SfSq y + f 2 5ql + k' 2 f 2 n5q 12 ) . 

By minimizing Eq. (|B5[) with respect to Sf, we find 
[2f 2 + {k/n) 2 } Sf = -2/0-70%, (B6) 
and solving for Sf 



Sf 



2/ 2 + (fc/^) 2 



fo 



go Sq y 
2/o 3 



(B7) 



where in the last step we kept only the leading and 
the next to leading order terms. Substituting back into 
Eq. (|B5l) gives 



S 2 T = 



dw 



1 - 3<Zo 2 ) Sql 



^ \ ?0c-2 K f 2 t -/2 

-) j2 S 1y + ^fo$% 



(B8) 



The first term in square brackets is the leading term. 
Neglecting the other terms, since q 2 , is monotonically de- 
creasing function of w the variation q y that minimizes 
the functional is a <5-function at the surface. Then the 
condition for the metastability is 



1-3<7 2 (0) = 0. 
Using Eq. (|A2j) we obtain 

cosh £ = Vd , sinh £ = \/b 
and substituting into Eq. (|A3|) 



(B9) 



(BIO) 



(Bll) 



To calculate the large-K correction, we expand the 
function qo(w) in the first term in square brackets in 
Eq. (|B8|) to linear order, while go and fo in the subleading 
terms can be simply evaluated at the surface. Setting 



£ ~ arccosh\/6 



1/4 



kK 3 /\ 



(B12) 



(B13) 



and using Eq. (|A3j) we find 



-C + w+lkA 5q 2 y + ^Sq' 2 



(B14) 



8 



The variational equation for Sq y derived from this func- 
tional has as solution the Airy function 



Sq y (w) = Ai 



5P 



1/3 



w — c + 



h 2 



Imposing the boundary condition 5q' y (0) 
that for a given k the lowest possible c is 

C = Zo [— ) + 4^' 

where 

z Q w 1.018793 

is the smallest number satisfying Ai'(— zq] 
minimizing c with respect to k we find 



4 \ 3/s /2 N 1 



and 



2 ^l/4 ^ 4 N 3/4 



(B15) 
= 0, we find 

(B16) 

(B17) 
= 0. Finally 

(B18) 
(B19) 



Substituting Eq. (|B12|) into Eq. ([A"5]l we obtain 
v/lO . 2c yiO . 0.3852 



6 3^3^ 6 
and from Eqs. ([BT3]) , (|BT7| . and (jBT8| 



(B20) 



160 \ 



1/8 



243 / 
0.9558k 3 / 4 . 



z 3/ V/ 4 



(B21) 



These results agree with those of Ref . [lOj . We compare 
these two formulas with numerics in Figs. [5] and [3J re- 
spectively. 

Finally, fixing the arbitrary normalization of the per- 
turbation by requiring Sq y (0) = 0, using Eqs. (IB15|) - 
(|B19[) . and restoring dimensions we find 



Sq y (x) = Ai 



10 

3 ° 



1/4 



A 



zo 



Ai[-z ], (B22) 



which shows that the "penetration depth" of the pertur- 
bation is of the order of the geometric average of coher- 
ence length and magnetic field penetration depth. This 
functional form is plotted in Fig.[5]for k = 50 along with 
the numerically calculated 5q y . 
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